function [flowData] = dicomAnalysis(dicomImg, inputAxes, instructionText)

iP = inputParser;

iP.addRequired('dicomImg', @isnumeric);

iP.parse(dicomImg);


% Perform color segmentation to eliminate everything except green ECG line
% and orange lines

RGB = im2double(dicomImg);
cform = makecform('srgb2lab', 'AdaptedWhitePoint', whitepoint('D65'));
I = applycform(RGB,cform);

% Define thresholds for channel 3 based on histogram settings
channel3Min = 19.089;
channel3Max = 70.537;

% Create mask based on chosen histogram thresholds
BW = (I(:,:,3) >= channel3Min ) & (I(:,:,3) <= channel3Max);

% Removes small artifacts
BW = bwareaopen(BW, 10);

% Perfoms image dilation to remove discontinuities (primarily in the ECG) to
% allow for better segmentation using region properties
[ dilatedBW ] = dilateEcg( BW );

% Isolates orange horizontal line
rp = regionprops(dilatedBW, 'PixelList', 'PixelIdxList', 'Eccentricity', 'Orientation');

% Removes ECG and orange axis information from the image
RGB(repmat(BW, [1 1 3])) = 0;

% Windows image to focus on velocity data
[ zeroAxis ] = segmentXAxis( rp );

[ xAxisStartPoint, xAxisEndPoint, xAxisBottom ] = findZeroAxisBoundsGUI( zeroAxis, dicomImg, inputAxes, instructionText );

[ windowedImg, leftSystoleBound, rightSystoleBound ] = getBoundingBoxGUI( dicomImg, RGB, xAxisBottom, inputAxes, instructionText );

%% Gets outline of max of velocity data

% Removes aliasing and returns velocity data
[largestVelConComp] = peakFinder (im2bw(windowedImg, 0.1), windowedImg);

% Identifies the coordinates of cells with value 1 in the velocity
% connected component
[ conCompVals ] = findBinaryCoord( largestVelConComp );

% Finds the maximum y value for each column of the velocity data (finds the
% outline of the data)
[ conCompMaxVals ] = findBinaryOutline( conCompVals );

%Gets datapoint range to be used for curve fitting
cla(inputAxes,'reset')
[ conCompMaxVals ] = getVelocityStartGUI( windowedImg, conCompMaxVals, inputAxes, instructionText );
[ conCompMaxVals ] = getVelocityEndGUI( windowedImg, conCompMaxVals, inputAxes, instructionText );

%% Smooth Data with Fourier Transform

%Fits curve to velocity data and smooths the fitted curve (ensures fit is
%good and prompts user for input)
cla(inputAxes,'reset')
flagFit = 0;
[ smoothedVelData ] = fitAndSmoothVelocityDataGUI( conCompMaxVals, flagFit, windowedImg, inputAxes, instructionText );

%% Find coordinates of R waves

% Select ECG
ecgSignal = bwareafilt(dilatedBW, 1, 'largest');

[ ecgPixelLoc ] = findBinaryCoord( ecgSignal );

% Removes duplicate Y (col) values to account for multiple pixels being in
% detected ECG line

[normEcgPixelLoc] = findBinaryOutline( ecgPixelLoc );
[unshiftedEcgPixelLoc] = normEcgPixelLoc; %saves a copy of the unshifted curve for later use with zero padding

% Ensures that ecg and x axis are same length for later indexing

[ normEcgPixelLoc ] = alignIndeces( normEcgPixelLoc, xAxisStartPoint, xAxisEndPoint );


%% Fits spline to ECG data to prevent discontinuities in data disturbing
% peak finding


[ fitResult ] = fitPoly( 'smoothingspline', normEcgPixelLoc );

smoothedEcgData(:, 1) = normEcgPixelLoc(:, 1);
smoothedEcgData(:, 2) = feval(fitResult, normEcgPixelLoc(:, 1));


%% Finds R wave peaks  
cla(inputAxes,'reset')
[ ecgPeakVals ] = getECGPeakValsManuallyGUI ( normEcgPixelLoc, unshiftedEcgPixelLoc, dicomImg, inputAxes, instructionText );

%Determines where to zero pad velocity data, and zero pads the velocity data
[ smoothedVelData, ecgPeakVals ] = addZeroPadding ( unshiftedEcgPixelLoc, normEcgPixelLoc, ecgPeakVals, rightSystoleBound, ...
    smoothedVelData, dicomImg );

%% Downsample data

% Number of points for downsampling
downThreshold = 128;

% Extracts ECG gated velocity, extracts axis information, smooths and
% downsamples data
[ downSampledVals, systoleVel ] = processBeatGUI( ecgPeakVals, ... 
    smoothedVelData, dicomImg, xAxisBottom, xAxisEndPoint, downThreshold, inputAxes, instructionText );

[ fitResult ] = fitPoly( 'smoothingspline', downSampledVals );

maxTimeVal = max(downSampledVals(:, 1));

refitVelData(:, 1) = min(downSampledVals(:, 1)):(maxTimeVal / downThreshold):maxTimeVal;
refitVelData(:, 2) = feval(fitResult, refitVelData(:, 1));

[ refitVelData ] = shiftFlowStartTimeToZero( refitVelData ); %shifts flow in time domain to start at time zero

% Sets negative values to zero to avoid issues with downstream processing
% for fourier analysis
refitVelData(refitVelData(:, 2) <= 0.1, 2) = 0;

%% Save relevant plots

[ ecgGatedVelPlot ] = plotEcgGatedVelOverRawData( windowedImg, systoleVel );

[ sampledEcgPlot ] = plotSampledEcg( BW, refitVelData, ecgPeakVals, xAxisStartPoint);
              

%% Create object for export

flowData.time = refitVelData(:, 1);
flowData.velocity = refitVelData(:, 2);
flowData.signalPeriod = length(systoleVel);
flowData.numberOfBeats = 1;
flowData.ecgGatedVelPlot = ecgGatedVelPlot;
flowData.sampledEcgPlot = sampledEcgPlot;
end
